{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 循环神经网络的变种\n",
    "在以上几节中已经完整地介绍了使用LSTM 结构的循环神经网络。这一节将再介绍循环神经网络的几个常用变种以及它们所解决的问题，同时也会给出如何使用TensorFlow来实现这些变种。\n",
    "\n",
    "### 8.3.1 双向循环神经网络和深层循环神经网络\n",
    "\n",
    "**1. 双向循环神经网络**\n",
    "\n",
    "在经典的循环神经网络中，状态的传输是从前往后单向的。然而，**在有些问题中，当前时刻的输出不仅和之前的状态有关系，也和之后的状态相关。这时就需要使用双向循环神经网络（bidirectional RNN）**来解决这类问题。例如预测一个语句中缺失的单词不仅需要根据前文来判断，也需要根据后面的内容，这时双向循环网络就可以发挥它的作用。双向循环神经网络是由两个独立的循环神经网络叠加在一起组成的。输出由这两个循环神经网络的输出拼接而成。下图展示了一个双向循环神经网络的结构图：\n",
    "<p align='center'>\n",
    "    <img src=images/图8.8.JPG>\n",
    "</p>\n",
    "\n",
    "从图8.8中可以看到，双向循环神经网络的主体结构就是两个单向循环神经网络的结合。**在每一个时刻t，输入会同时提供给这两个方向相反的循环神经网络。两个网络独立进行计算，各自产生该时刻的新状态和输出，而双向循环网络的最终输出是这两个单向循环神经网络的输出的简单拼接。**两个循环神经网络除方向不同以外，其余结构完全对称。每一层网络中的循环体可以自由选用任意结构，如前面介绍过的简单RNN、LSTM均可作为双向循环网络的循环体。双向循环神经网络的前向传播过程和单向的循环神经网络十分类似，这里不再赘述。更多关于双向神经网络的介绍可以参考Mike Schuster和Kuldip K. Paliwal发表的论文[Bidirectional recurrent neural networks](https://maxwell.ict.griffith.edu.au/spl/publications/papers/ieeesp97_schuster.pdf)。\n",
    "\n",
    "**2. 深层循环神经网络**\n",
    "\n",
    "深层循环神经网络（Deep RNN）是RNN的另外一种变种。**为了增强模型的表达能力，可以在网络中设置多个循环层，将每层循环网络的输出传给下一层进行处理**。在图8.2描述的单层循环网络中，每一时刻的输入$x_t$到输出$o_t$之间只有一个全连接层，因此在$x_t$到$o_t$的路径上是一个很浅的神经网络，从输入中提取抽象信息的能力将受到限制。图8.9给出了深层循环神经网络的结构示意图，可以看到，在一个L层的深层循环网络中，每一时刻的输入$x_t$到输出$o_t$之间有L个循环体，网络因此可以从输入中抽取更加高层的信息。**和CNN类似，每一层的循环体中参数是一致的，而不同层中的参数可以不同。**\n",
    "<p align='center'>\n",
    "    <img src=images/图8.9.JPG>\n",
    "</p>\n",
    "\n",
    "为了更好地支持深层循环神经网络，TensorFlow中提供了MultiRNNCell类来实现深层循环神经网络的前向传播过程。以下代码展示如何使用这个类："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "# 定义一个基本的LSTM结构作为循环体的基础结构。Deep RNN也支持使用其他的循环体结构。\n",
    "lstm_cell= tf.nn.rnn_cell.BasicLSTMCell\n",
    "\n",
    "# 通过MultiRNNCell类实现深层循环神经网络中每一个时刻的前向传播过程。\n",
    "# 其中number_of_layers表示有多少层，也就是图8.9中从xt到ht要经过多少个LSTM结构。\n",
    "# 注意从TensorFlow1.1版本起，不能使用[lstm_cell(lstm_size)] * N 的形式来\n",
    "# 初始化MultiRNNCell，否则TensorFlow会在每一层之间共享参致。\n",
    "stacked_lstm = tf.nn.rnn_cell.MultiRNNCell(\n",
    "                [lstm_cell(lstm_size) for _ in range(number_of_layers)])\n",
    "                \n",
    "# 和经典的循环神经网络一样，可以通过zero_state函数来获取初始状态。\n",
    "state = stacked_lstm.zero_state(batch size, tf.float32)\n",
    "\n",
    "# 和8.2节中给出的代码一样，计算每一时刻的前向传播结果。\n",
    "for i in range(len(num_steps)) :\n",
    "    if i > 0: tf.get_variable_scope().reuse_variables()\n",
    "    stacked_lstm_output, state = stacked_lstm(current_input, state)\n",
    "    final_output = fully_connected(stacked_lstm_output)\n",
    "    loss += calc_loss(final_output, expected_output)\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从以上代码可以看到，在TensorFlow中只需要在`BasicLSTMCell`的基础上再封装一层`MultiRNNCell`就可以非常容易地实现深层循环神经网络了。\n",
    "\n",
    "### 8.3.2 循环神经网络的dropout\n",
    "6.4节介绍过在CNN上使用dropout的方法，通过dropout可以让CNN更加健壮（robust）。类似的，在RNN中使用dropout也有同样的功能。而且，**类似CNN只在最后的全连接层中使用dropout，RNN一般只在不同层循环体结构之间使用dropout，而不在同一层的循环体结构之间使用。**也就是说从时刻t-1传递到时刻t时，RNN不会进行状态的dropout；而在同一个时刻t中，不同层循环体之间会使用dropout。\n",
    "\n",
    "图8.9展示了RNN使用dropout的方法。假设要从t-2时刻的输入$x_{t-2}$传递到t+1时刻的输出$y_{t+1}$，那么$x_{t-2}$将首先传入第一层循环体结构，这个过程会使用dropout。但是从t-2时刻的第一层循环体结构传递到第一层的t-1、t、t+1时刻不会使用dropout。在t+1时刻的第一层循环体结构传递到同一时刻内更高层的循环体结构时，会再次使用dropout。\n",
    "<p align='center'>\n",
    "    <img src=images/图8.10.JPG>\n",
    "</p>\n",
    "\n",
    "在Tensorflow中，使用`tf.nn.rnn_cell.DropoutWrapper`类可以很容易实现dropout功能。以下代码展示了如何在TensoFlow中实现带dropout的RNN:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "# 定义LSTM结构。\n",
    "lstm_cell = tf.nn.rnn_cell.BasicLSTMCell\n",
    "\n",
    "# 使用DropoutWrapper类来实现dropout 功能。该类通过两个参数来控制dropout的概率，\n",
    "#   一个参数为input_keep_prob，它可以用来控制输入的dropout概率；\n",
    "#   另一参数为output_keep_prob，它可以用来控制输出的dropout概率。\n",
    "# 在使用了DropoutWrapper的基础上定义MultiRNNCell。\n",
    "stacked_lstm = tf.nn.rnn_cell.MultiRNNCell(\n",
    "              [tf.nn.rnn_cell.DropoutWrapper(lstm_cell(lstm_size))\n",
    "              for _ in range(number_of_layers)])\n",
    "\n",
    "# 和8.3.1节中深层循环网络样例程序类似，运行前向传播过程。\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.4 循环神经网络样例程序\n",
    "在以上几节中已经介绍了不同RNN的网络结构，并给出了具体的TensorFlow程序来实现这些RNN的前向传播过程。这一节将以时序预测为例，利用RNN实现对函数$sin(x)$取值的预测。图8.11给出了$sin$函数的函数图像。\n",
    "<p align='center'>\n",
    "    <img src=images/图8.11.JPG>\n",
    "</p>\n",
    "\n",
    "下面的篇幅将给出具体的TensorFlow程序来实现预测正弦函数$sin$。因为RNN模型预测的是离散时刻的取值，所以在程序中需要将连续的$sin$函数曲线离散化。所谓离散化就是在一个给定的区间[0, MAX]内，通过有限个采样点模拟一个连续的曲线。比如在以下程序中每隔SAMPLE_ITERVAL对$sin$函数进行一次来样，采样得到的序列就是$sin$函数离散化之后的结果。以下程序为预测离散化之后的$sin$函数:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train step: 0, loss:  0.4196214\n",
      "train step: 1000, loss:  0.0010739247\n",
      "train step: 2000, loss:  0.0002879769\n",
      "train step: 3000, loss:  5.4690085e-05\n",
      "train step: 4000, loss:  1.3475828e-05\n",
      "train step: 5000, loss:  4.33033e-06\n",
      "train step: 6000, loss:  5.110061e-06\n",
      "train step: 7000, loss:  4.500361e-06\n",
      "train step: 8000, loss:  3.534146e-06\n",
      "train step: 9000, loss:  3.5462117e-06\n",
      "Root Mean Square Error is: 0.002050\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xl8VOXd///XJ3tCErISErISEsjKFnaUHQFRXNCKda1b7W2r7W3v2l/bu613/dba3mr1tlUsrrVuKIqIIqvKJgQFQhJCVpIQlmyE7Ntcvz8y2ICBJMxkzszkej4e88jMmevMeU8mySfnOudclyil0DRN07SzXIwOoGmaptkXXRg0TdO0c+jCoGmapp1DFwZN0zTtHLowaJqmaefQhUHTNE07hy4MmqZp2jl0YdA0TdPOoQuDpmmadg43owNcipCQEBUbG2t0DE3TNIeyb9++KqVUaG/tHLIwxMbGkpmZaXQMTdM0hyIiR/vSTnclaZqmaefQhUHTNE07hy4MmqZp2jl0YdA0TdPOoQuDpmmadg6rFAYReUlETonIoQs8LyLyjIgUiMhBEZnQ7bnbRSTffLvdGnk0TdO0S2etPYZXgEUXeX4xkGC+3Qv8HUBEgoDfAlOAycBvRSTQSpk0TdO0S2CV6xiUUl+ISOxFmiwDXlNd84juFpEAEQkHZgMblVI1ACKyka4C86Y1cg0qJhOczKK1bD9NNccwiRtDgiPwipkMIQkgYnRCTRs41YV0lH9Nw8kiOjs68AgYjm9UGjJiIri4Gp3O4djqArcRQFm3x+XmZRda/h0ici9dextER0cPTEpH1HCKqi3P4nXwn/h21OAJeJ7X5JRnNO3ptzBiwQPgMcSIlJpmfR2t1O96mbbdLxLcWIAbEHBek3rXAKrilhG99Be4BvT4p0Xrga0KQ0//rqqLLP/uQqVWAisBMjIyemwzqJg6qd7yDD47niDI1MxWNZ7C0LvxjJuO/7AY3KSTxlPFdJTsZMypT8jY+/+o3fd36mf/D9GX32p0ek2zSFP2J7R9+CABbSfZbxrJ+3734Ro3k6EjRuPh4UFL7TFaj+4lrPwz5uS/TufTb1CU/CNGXfdrxO38f52089mqMJQDUd0eRwIV5uWzz1u+zUaZHJaqP8Hxl24honYvX6rxFE38Fcvmz2Kej8d5LROAhTS0/pqPN64lJvMxUrc8QM7BdSTevQo3L18j4mvapetsp+LtnxJx5HWOmUbw7sj/Zf6VN3FP6Pk/y3HATDpND7F9byZq0++ZnfMMJQWfEviDtxk6fKQR6R2GdHX7W+GFuo4xrFNKpfbw3JXAA8ASug40P6OUmmw++LwPOHuW0tfAxLPHHC4kIyNDDdaxklpP5NH4j6vwbq/lrdCfsOTWhwkb6t2ndesam9n9yq9YcOolij1GEXrvB/iHRg5wYk2zDtVcy7GVNxJZu4f3PJaRcPOfSY8N69O6HZ0mtn7wElMP/poOcefMdf8iJv2yAU5sf0Rkn1Iqo9d21igMIvImXf/5hwAn6TrTyB1AKfW8iAjwf3QdWG4C7lRKZZrX/QHw/5lf6jGl1Mu9bW+wFoam0v20v3I17Z2KzRP/xo1XLUUu4aDyjo9fY/yeh6lyDcX7nk8IDdfHbDT7pppPc+zZxQxrzOOt4Q+z/K7/wsej/x0euVn78H/vJvypp+LKfzJ60vwBSGu/bFoYbG0wFobGE/m0r5xPU6creQtfZ86MGRa93oHtH5Ow8U5OuoYRcP9nBIaGWympplmXam2g/K8LCWs8zJqEP3Lj9++9pH+IzjpWko/p1asJMtVQfs1qRo8fPHsOfS0M+spnB9Bx5hT1/7ga1dnBESsUBYCxM6/k6BUvEdF5nGMrb6CpuckKSTXNykwmClfeQkRjDmsTHrO4KACMiE3A+55PaHDxJejDWyktPmKlsM5DFwY7pzrbKV35PQLaK8mc/ndmW6EonJU0fSl5U/4fqe1Z7Pvb3SiTyWqvrWnWcOStXzCqeisfDX+A67//Q4uLwlkhEbG0fe8tvGil9fUbqDtzxiqv6yx0YbBzuf96hJENX7Np5CMsuOJqq79++pJ72R/zAy6r/5ht7/zV6q+vaZeqbPcaEo+sZJP3Ypbc/XurFYWzosZM4sSC50gwlbD/xfswmRyvW32g6MJgx45lriO58B9sGbKEJbf+54BtZ+xtf+aI9zgm5/6RA/v3Dth2NK2vmmqO4bvhQfKJZtx9L+DhPjBXLyfMuI6DsT9gVv16trz73IBswxHpwmCnWupr8fj4QYqIJPWuv+PiMnBDWoirGyN+8Bod4oHnh/fR0NQ8YNvStF4pRdlLd+BtaqJ+6UpCAoYO6ObSbn2CQq9UJuc8Rl5+3oBuy1HowmCn8l59gEBTDTUL/8qwoPMv9Le+IaExVM15nDGqkC9f+92Ab0/TLiRn/XOMbtjDlyMfYkLGtAHfnri6E3LLKjykg9q3H6C1vWPAt2nvdGGwQyW7P2Rs1Tq+DPs+GdNtd551/KxbOBI4iznHV5H59T6bbVfTzmqsOU7k3j+S5ZrCrJsfsdl2h0aOoXTcz5jasYdN7/zNZtu1V7ow2JmO1mY8PnuEo0Qw4dbHbb796Fueo0PckY8f0v85aTZX+PqP8VLNuC77Kx7uthqxp0vi1f9FiVcS0448QemxCptu297owmBnvnnnMSJMFZyY8ShD/Ww/lpFXcBTHMv6LiZ0H2brmRZtvXxu8CvesJ712IzvDbyc5fZLtA7i44r/8WQJoIOetX+OIF/9aiy4MduT08WJSCleS6T2DyfOXG5Zj9JKfUOYxktTs/6W88qLDVmmaVShTJy6f/YrjhDDx+783LEfQqEkcGXEt8858wM7dOw3LYTRdGOxIyTuP4KJMBF/3F6ufs90vLq54Xfk4kVLJN+/avjtLG3wOrl9JXEcRxWMfxs/P39As8Tf9iVbxwn3Tr2jvHJwXferCYCdKczNJr9nA3rAbiEtINjoOoWOvID9gJrNPvsaRoiKj42hOrLW5nvB9fybPNYEpV99rdBzc/YdxfNyPmdz5Dds+XW10HEPowmAnaj76DQ14k3zDb42O8q3h1/8Zb2mlaM0fjI6iObGD7/6RYaqalrmP4upqH9NwjrryIapcQgjL/DNNre1Gx7E5XRjsQMG+TYxr2smh2DsIDh1udJxv+UUlUzD8SmafWcu+Q7lGx9GcUGNdDaOLXuFr76mMnbHE6DjfEndvGqf+jHR1hC1rXzc6js3pwmAHOj77PVUMJe36Xxgd5Tvirvs9btLJqU/0sQbN+rLXPIE/jfgs/LXRUb4jZt69nHKLYFT209Q3txodx6asUhhEZJGI5IlIgYh856oUEXlKRPabb0dE5HS35zq7PbfWGnkcyZE9GxjTepDc+Lvx8x/4K5z7y3NYPAURy5jbsI4D2dlGx9GcSENdDaNLXuNr72mMscc5EVzdaZ35C8ZwlJ0fvWR0GpuyuDCIiCvwHLAYSAZWiMg5R0+VUj9VSo1TSo0DngXe7/Z089nnlFLWHz7UzrVs+TM1+DP+2oeMjnJBMdf8Ny4CVZ/80egomhPJXvMEQ2nEZ8H/13tjg0RdfisVblHE5DxPc+vgueDTGnsMk4ECpVSRUqoNeAtYdpH2K4A3rbBdh3fkmy9Jb9nLkbjb8fU19hS9i/EeNpIjw5cys/5Tco8UGB1HcwKNZ2oZU/Ia33hPZcyEy42Oc2EurrRO+TFjKOGLT982Oo3NWKMwjADKuj0uNy/7DhGJAeKALd0We4lIpojsFpFrrJDHYTRufoIzyoeUZT81OkqvYq76Be7SQdEnTxkdRXMCh9b9H0NpxGue/R1XO1/cnDupcgkhdP9ztHZ0Gh3HJqxRGHq6EutC15LfBKxWSnX/7kab5yC9GXhaROJ73IjIveYCkllZWWlZYjtQkX+AsfVfkjXie/gFBBsdp1e+I5IpDJrFjJo1lFScMjqO5sA62tuIzn+VHPdUkjLmGh2nd24e1I27jwkqm+1b1xudxiasURjKgahujyOBC41AdRPndSMppSrMX4uAbcD4nlZUSq1USmUopTJCQ0MtzWy4ik//lzbcGHXVw0ZH6bPgK35OgDSS+/H/GR1Fc2AHN75OuKqkZdL9Rkfps5FX3M8Z8cNnz7ODYgwlaxSGvUCCiMSJiAddf/y/c3aRiIwGAoFd3ZYFioin+X4IMAPIsUImu3a66gSpVev5JnARYeGRRsfps6DRMyn0SWds+RvUNTQZHUdzQMpkwvfr5ymTCMbOW2F0nD4TTz/K4m9mStse9u//2ug4A87iwqCU6gAeADYAucA7SqlsEXlURLqfZbQCeEudW26TgEwROQBsBR5XSjl9YchZ9yxe0k7Ywp8YHaXf3Gb+hAipYs+nrxkdRXNAh/dsILHjCMeS7rSbq5z7Kn7Jg3SKC1VbnX8KUHHE3aKMjAyVmZlpdIxL0tbaSu0fx1DlGUPKL7cZHaf/TJ2cfCyZClMQab/egZurvkZS67v9TywmpikLr5/n4j3Ez+g4/Zb7fzcyovILau87QExEmNFx+k1E9pmP6V6U/q22sYObXieMGjqn/NDoKJfGxZXTKbcxXuWwa9cXRqfRHMjxo3mkN+7i8IjlDlkUAMIWPIi/NHPokxeMjjKgdGGwMf/9L1Iu4aTOusHoKJds1BX304IH7bueNzqK5kCOfvYcCohb9IDRUS5Z0OgZlHglkVz2Js1OPLieLgw2VJy1ncT2wxwddQsuDta/2p3rkCAKhy9hWsNmSsuPGR1HcwBtrS0kHPuAg0OmMjxqlNFxLNKRcQ9xVLB383tGRxkwujDYUNW252lWHqQsus/oKBYbvuAneEsb+Rv+bnQUzQFkbXqdYOpwmXS30VEsFj/7FmokAK/9zjt+ki4MNlJfV0NK1WdkBcwjINjxr8MIjp9IvlcaiWXv0NruvLvUmnX4HHiFYxJG2uXXGh3FYuLmSWn0tUxs3UN+QZ7RcQaELgw2cmjDS/hIKwGXGz9DlbV0TLiTKE6y7/OPjI6i2bGjOXtJajtEadxNDt2F2l3cgvtxFcXRTS8aHWVA6MJgA0opgvPepNg1lsQJc4yOYzWjZ62gniGofa8aHUWzYye3/p1W5c7oRQ56Jl4PhkaO5ojPBMYc/4Cm1jaj41idLgw2UHhgO4mdBVQmrgDpaWgpx+Ti6UNxxFIymnZwtLzc6DiaHWppqiepcj0Hh84maFiE0XGsyiXjTiKlkr1b3u+9sYPRhcEGar98kWblweiFjn/g7Xwj5t2Hp7RTsMl5D8Rply578xv40Yz35DuMjmJ18ZfdSB1+eB5wvqk/dWEYYC0Np0mu3sDBgLkMDQwxOo7VBcdPpMgjkeij72HqNBkdR7Mznofe4piEkTxtsdFRrE7cvSiOvJqJzbsoLztqdByr0oVhgOVueYMhtOA95U6jowyY+uQVJKgSsjK3GR1FsyMnS4+Q3LKfo5HLnOag8/ki5t6Hu3RS7GQHoXVhGGCe2W9TJuGkTllodJQBM3renTTjQePul42OotmRks3/wEUUMXOdrwv1rGEjx5LnkUJU6Rqn2mPWhWEAHT+aR3LrAcqirsbFiQeb8/ILJCdwLuk1G2loqDc6jmYHlKmTqNIPyPIYx4i40UbHGVCNSTcQq8rJ2fe50VGsxnn/WtmB4s1dB2Tj5t1lcJKB5zf5VnylmawtejpvDfL2bCBCnaQ55Sajowy4MfNuo1W5c3q38xyE1oVhgCiTiaiyteR4pBMe49z/MQEkTFnMKQnGM/tdo6NodqBh96s0KG9S5n3f6CgDzsc/mMNDZ5Jc/RnNzS1Gx7EKXRgGyJGvtxClKmhIctxRVPtDXFwpG7GU9JZMKo6VGh1HM1BTfS3JtVvJDpzHEF9/o+PYhGfGzQRJPQe2rTY6ilVYpTCIyCIRyRORAhF5pIfn7xCRShHZb77d3e2520Uk33y73Rp57MGZ3a/TrDwYM/cWo6PYTMSsO3ETE0Vb9exug1nuln/hI634Tr3N6Cg2kzjtGmrxxzXrLaOjWIXFhUFEXIHngMVAMrBCRJJ7aPq2Umqc+fYP87pBwG+BKcBk4LciEmhpJqO1tzYxumojWf6X4z80yOg4NhOeMJ4it5EMK/7A6CiagTxz36OCYSRNWmB0FJtxcfegYPhi0ht3c7r6lNFxLGaNPYbJQIFSqkgp1Qa8BSzr47pXABuVUjVKqVpgI7DICpkMlffFu/jTiMv4m42OYnNVI68lsTOfksPfGB1FM0DtqXKSmr+mOGKxU5+J15Og6bfhKe3kbXH8g9DW+ORGAGXdHpebl53vehE5KCKrRSSqn+s6FNOBdzhFIOkzrzY6is2NnHM7nUo4sV0PrDcYFWx9HVdRDJ9+q9FRbG5k6nRKXKIYmu/4xxmsURh6GhVOnff4IyBWKZUObALO/tXoy7pdDUXuFZFMEcmsrKy85LADrbGuhjH1uzkSsgAPD3ej49hcSHgM2V4TiDm2DmXqNDqOZmNDCz6g0CWWkSm9zjfvdMTFhfKoqxnTlkNlqWPP02CNwlAORHV7HAlUdG+glKpWSrWaH74ITOzrut1eY6VSKkMplREaar8T3Rze9i88pIPAySuMjmKYpjHLCVeVFGZuMjqKZkMnjh4msf0wJ6KXIk40inB/RF/etadU8sUbBiexjDUKw14gQUTiRMQDuAlY272BiIR3e3g1kGu+vwFYKCKB5oPOC83LHJZ77hqOSRhJE2cbHcUwSXNX0Kw8qN3rHGdoaH1Tuq2rIyBm9uA5G+l80fFJHHYbTXCJY09eZXFhUEp1AA/Q9Qc9F3hHKZUtIo+KyNlO9p+ISLaIHAB+AtxhXrcG+B+6iste4FHzModUdfIYyc1fUxq+aNAdeOtu6NBAsn2nEl+5mc4OPe3noKAUw49+RLZ7CpGxzn9B58VUx17FyI4ijuXvNzrKJbPKXy+l1HqlVKJSKl4p9Zh52X8rpdaa7/9SKZWilBqrlJqjlDrcbd2XlFKjzDeHHoWt4PM3cBMT4TOc/2rPXqVcRxB1HN79idFJNBsozd1DtKmMuvhrjI5iuJGzb8GkhPLtjtudNHj/rR0A/gVrOeoSRWzyZKOjGC5l1nIalRdNX79jdBTNBk5uf5125UrCHP1PUXhkHLkeqYSXrQfV47k0dk8XBiupOlbEmNZDHItc4lTTd14q7yF+5PjPIKFmK53tzjcnrvZvymQi6vinZHtPJDTM4c82t4r6hKuJNpVTkrPH6CiXRBcGKyn+4g1cRBGhu5H+LeV6Amggf7djH4jTLq7gm20MV5W0jtbdSGfFz7qZDuVC5S7HHG1YFwYrCSj6iHyXeGJHjzU6it1Iufwazigfmr/RI646s5q979KmXBkz2/mH2O6r0LBIDnmNJ7LCMbuTdGGwgqrSwyS053Ei+kqjo9gVH58hZPtfxqiazzG1OcdwxNq5lMlE9MmN5PhkMDQw2Og4dqV+1NWEm05Snr3d6Cj9pguDFRz94p8ARF02+MZG6o2kXocfTRTu/tDoKNoAKDrwJeGqktbEwTf8S28SLr+JVuVG1a5/GR2l33RhsIKhRzeQ55pAbHyS0VHsTuply6hVvrTo7iSnVLXnHdqUK6Mvv9HoKHZneNhwDnpOIOL4RofrTtKFwUKnyosY1X6E6qiFRkexS74+3mT5zyK+9ktMrU1Gx9GsSJlMRJ/4jFzvCQQEDzM6jl1qGnUlw0yVnMjZaXSUftGFwUKF27vO0x8xdXDM1HYpXFOvxocWiveuMzqKZkXFWTsIV6doTrjK6Ch2K37mjbQrV0589bbRUfpFFwYL+RR9SpnLCGLGjDc6it1Km3kVdWoIjd+8b3QUzYqqvnqbduXK6Fn6bKQLiYyIIMtjLGHlnzlUd5IuDBY4cfI4ya0HORUx3+gods1/yBAO+U4ntvoLVIe+2M0ZKJOJyOMbyfEaT2BImNFx7NqZuMWEm45zsmCf0VH6TBcGC+R98R7u0knYlOVGR7F/yVfjTyPFmZ8anUSzgpLsXUSoEzQlLDU6it0bOfN7dCqhYqfjdCfpwmABj4L1VEsQkSkzjY5i91JnXkOj8qTu6/eMjqJZQeXut+lQLiRcrruRehMdHUO2ewohpY7zT5EuDJfoVHUtY1syqQibAy7629iboUP9OTRkKjGntqI6O4yOo1lCKUYc/4xsr3GEDAvvvb1GbcxiojpLqSo+aHSUPtF/0S5Rzo61+EgrgRnXGR3FYbQlLiWIOsoObDU6imaBozl7GGE6TlO87kbqq+gZ3wOgbIdjTF6lC8Mlkrz1NODDiHELjI7iMMZcdh2typ2aTMefLH0wO/7VakxKGKUvauuzuJEJZLuMJvCoY3QnWaUwiMgiEckTkQIReaSH538mIjkiclBENotITLfnOkVkv/m29vx17dGZpmbSGnZQEjwTcfM0Oo7DCA0O4YDnREac2ORQp+5p5wo5tok8j2RCh0f13lj71snIK4htL+RMRb7RUXplcWEQEVfgOWAxkAysEJHk85p9A2QopdKB1cAT3Z5rVkqNM98cYsCVg7s+I0jq8U5ziLh2pSF+CaGmKk7lOdaVoFqX0sJcRnUWUR97hdFRHM5w80Wwpdvtfyhua+wxTAYKlFJFSqk24C1gWfcGSqmtSqmz4yHsBiKtsF3DtGWtpRV3Yqcs672xdo74GctpV66c3K1ndnNEZbu6ugFjZugr/ftrzJg0DstIhhSuNzpKr6xRGEYAZd0el5uXXchdQPeJgL1EJFNEdovIBWf6EJF7ze0yKysrLUtsgdb2DhJrP6fYLwNXb3/DcjiqmMgRHHBPJ6TMsa4E1br4H/2MEtcYwmLP7xTQeuPiIpQOX0Bcay4t1UeNjnNR1igMPc1j2eNvvIjcAmQAf+62OFoplQHcDDwtIvE9rauUWqmUylBKZYSGhlqa+ZId3LedSKmEJD0+zKWqiV5EeGcFZ0r2Gx1F64fKkxUktR2iMmKe0VEcVuDE6wEo2W7fow1bozCUA92PQkUCFec3EpH5wK+Aq5VSrWeXK6UqzF+LgG2AXQ86VPf1B5iUEDfjeqOjOKyoqcvpVELZTsc4dU/rUrDjfdzEROgk/bN/qcaOm0QRI3A58knvjQ1kjcKwF0gQkTgR8QBuAs45u0hExgMv0FUUTnVbHiginub7IcAMIMcKmQZEp0kRfWoLRT6peA4dbnQchzUmIZ4slzH4H91odBStH9zz13NKgolJnW50FIfl4eZCYdAs4hq/oaOhxug4F2RxYVBKdQAPABuAXOAdpVS2iDwqImdP2/kz4Au8e95pqUlApogcALYCjyul7LYw5OQcIJGjtMQvMTqKQxMRTkbMI6qtkOZTRUbH0fqgvv4MyU2ZlIbOQfSV/hbxSV+GO52U7F5jdJQLssonrJRar5RKVErFK6UeMy/7b6XUWvP9+UqpsPNPS1VK7VRKpSmlxpq/rrJGnoFyYk/XsNExM/UZGZYKzbgWgJId9t3XqnXJ3fERPtKK3zh9Jp6lxk6dy0kVSFv2R0ZHuSBd+vtIKcWw8o2Uuo/Eb3iC0XEcXnr6BAqIwi3fvvtatS6dueuox4dRkxYZHcXh+Xp5kOM3g9jaXaj2ZqPj9EgXhj4qLCkhtTOX0zF6CAxrcHN1oThkNnFNB2ivrzI6jnYRrW1tjD79JQVDZ+Dq7mF0HKcgSVfiQwtl+zYYHaVHujD00dGdq3EVReRUPT6MtfiOvRo3TBTv0jO72bOcrzYSJPW4p+hTtK0lZfpS6pU3dd/Y53EGXRj6yK/kM065DCMofqLRUZzGuMlzOKGCaM/Wc0Hbs4aDH9Km3Bg1XR9fsJbQQH8Oek8i8tQ2MHUaHec7dGHog4pTlYxt+4YTEfNBerqeT7sU3p7uHB46k5F1u1FtTb2voNmcqdNEbOXnHBkyES/fAKPjOJXWkYsIVKc5mbvd6CjfoQtDH+RtX4OntBOi516wOtekpXjTSvFefRDaHuVm7SGKE3QkLjY6itMZNfM62pQrJ/fY36yGujD0gXv+J9SJHxFpc4yO4nRSZ1xJvfKm4cAHRkfRelCV+R4mJcTrU7StLjoinCz3dELK7W8Yel0YelF7poH0pt2UhswCVzej4zidQH9fsnymEFX5uV32tQ5mSimGV2ym0DMJvxCHHhDZbtVGLyCi8xinSw8ZHeUcujD04tDO9fhLE75jLzjwq2ahtlGLCVR1HDv0hdFRtG6OFh1htKmQM7ELjY7itCKnLgegdJd9XeipC0MvOnM+ohlPYidfaXQUp5Vo7mut3Gt/fa2DWal57oXo6foU7YEyOiGRHBmFb7F9Xc+gC8NFNLe2k1T3JUVDpyIePkbHcVoRYcPI8hjLsIrNdtfXOpj5H/2MMtcoQmNTjI7itESE4+HzGNl6mKaqUqPjfEsXhos4+NUWwqQW9xQ9hedAOxO9kIjOCqpLDhodRQNOnDxOSlsWlSP03AsDLWhi17hhxTtWG5zk33RhuIiGAx/SgQtx0/VpqgMtelrXGP9ldtbXOlgVbH8fd+nUcy/YQNrYyRwlHFc7mqNBF4YL6Og0EVe9jSKf8bj7Bhkdx+mNjE8gxyUBvxI9R4M98Cj4hCoJJCplptFRnJ6bmytFwbOIb9xHe2Ot0XEAXRguKOtAJiM5RkeinnvBFkSEk+FziW87TEOl/fS1DkZ1Z+pJadpDaegc0HMv2IR32tW400nRLvu4nscqn7qILBKRPBEpEJFHenjeU0TeNj//lYjEdnvul+bleSJyhTXyWEN1ZtcZMnH6wh6bCTF3WxTZUV/rYJSz8yOGSCv+4/SxNVtJnzqfKjWUdjuZo8HiwiAirsBzwGIgGVghIsnnNbsLqFVKjQKeAv5kXjeZrqlAU4BFwN/Mr2copRThxzdT7DEa75AYo+MMGslpkyhlOG5H1hsdZVDrzFlHA96MzNDDYNiKj5cnOX4ziKvdiWpvMTqOVfYYJgMFSqkipVQb8BZw/jCMy4BXzfdXA/NERMzL31JKtSqlioEC8+sZKi//CCkqn4Y4u9mBGRRczXM0JDR+TWuDffS1DjYtrW2MqdtOUcB0XDy8jI4zqEjSEobQzNF9nxkdxSqFYQRQ1u1xuXlZj23Mc0TXAcF9XNc6K/xpAAAgAElEQVTmjpkv7InSF/bYnO/YZbhLJwU77aOvdbA5tGczIVKHW7Kee8HWkqdfRaPy5PQ3xv/sW6Mw9DQO9flXKV2oTV/W7XoBkXtFJFNEMisrK/sZsX8Cyj6jwnUEAdGpA7od7btSJs+nWvnTkWMffa2DTcOBD2lXroyaroeAsbXgwAAOeWcQeWormEyGZrFGYSgHoro9jgQqLtRGRNyAoUBNH9cFQCm1UimVoZTKCA0NtULsnpUdO0Z6exaVkQv13AsG8PL04PDQmcSf3ompzfi+1sGks9NEXNU2CoZMwMM30Og4g1LzyEWEqBqOH95laA5rFIa9QIKIxImIB10Hk9ee12YtcLv5/nJgi1JKmZffZD5rKQ5IAPZYIdMly9+xBnfpZPhkfWGPUVySluJLM4WZnxodZVDJycokhuN67gUDjZpxPR3KxfA5GiwuDOZjBg8AG4Bc4B2lVLaIPCoiZ893WwUEi0gB8DPgEfO62cA7QA7wKfAfSilDx172LvyEagkiLGmGkTEGteQZXX2t9fuN72sdTM4OYhg3Q5+ibZTIESM45J5KUPkmQ3NYZYIBpdR6YP15y/672/0WoMefNqXUY8Bj1shhqcqa06S37KUwfCnB+sIewwz192ePzyTiKj9HmToRF8PPYHZ6SimGH99MoccY4kOjjY4zqNVGzWdc8ZPUlOUSFJVkSAb916+b3B1dF/YMHa8PvBmtbdRiQlUN5dk7jI4yKBQWHiHZlE99rD5F22gRU7q6sUt3Gnehpy4M3UjeOhrwIXqi/uUw2qiZ19GhXDi1532jowwKZbu6vs9R05cbnERLHJ1CvsTiU2TcMTZdGMwamltJqd9BSdAMxM3T6DiD3vCwCLI90hhWYWxf62DhX9p1inZwTJrRUQY9EeFY2FziW7JprDluSAZdGMyydm0gSOrxTD3/om3NKHUxC4nqLONUcZbRUZxaxcmTpLUd4FTEPH2Ktp0InHgNrqIo2mHM2Um6MJi1ZK2lFXdGTtOFwV5ET+vq1ijbpQfVG0hHvnwPD+lkmD5F226kjJ9JBSFInjHjhunCALS1d5JQ+zlFfhm4evsbHUczixk5miMu8fiXGD92jDPzLPiEGgkgIuVyo6NoZm5urhQGzWJUwx7am+ttvn1dGICsr3cQySkYfaXRUbRuRIQTEfOIb83lzKlyo+M4pdq6etKa91Cm516wO15pV+FFO4W7bT88jP5JAE5/swaTEuJm6DMy7E1oxrW4iKJQz9EwILJ3rMNXWvAfr7tQ7U3q1MXUqSG0ZOnCYHMmkyLq5BaKvFPxCgw3Oo52ntFpUznGMNzz9RwNA8F0+CMa8SJWz71gd7y9vcjxm0ZszZeoznabbnvQF4ac3CwSVQkt8YuMjqL1wMXVheKQOSQ27qOl4bTRcZxKc2s7Sea5F8Rdz71gj1TiEgKop/jrLTbd7qAvDCe+6uqiiJmh516wV0PGXo2HdJCv52iwqoO7NxIqdbin6LkX7NWYy66lVblT8/Uam253UBcGpRSh5ZsodY/DLyLR6DjaBSRPWUCN8qMje53RUZxK08G1tONK/PRrjY6iXUBQYBCHvMYz4sQWUD1OVTMgBnVhKCopIbUzh9PRC42Ool2Ep4cneUNnMqpuB53tbUbHcQodHZ3EV2+lYMhE3IfouRfsWfPIKwhXJzmWt89m2xzUhaFk52pcRRE5TXcj2TvX5KX40UTBXj1HgzVkH9xLNCfoTFxidBStF3HTl2NSwvGvbHdm3qAuDP4ln3LCZThB8RONjqL1ImnGVTQrDz1Hg5VUZ3YNmjdSn6Jt90ZExXLYbTSBZbYbN2zQFoZjJ06S3rafkyPm6/FhHICf31ByfCYRXbkNZfB8uI6u+9wLPiFRva+gGa46cgHxHfnUVBTZZHsWFQYRCRKRjSKSb/76nc5KERknIrtEJFtEDorI97o994qIFIvIfvNtnCV5+uPI9vfxlA5CJ+v/mBxF66hFDFPVHD2k52iwxJH8PJJVAQ1xenh5RzF8ynUAFO941ybbs3SP4RFgs1IqAdhsfny+JuA2pVQKsAh4WkQCuj3/c6XUOPNtv4V5+syr4GPz+DCzbLVJzUIJly2nUwmn9uo5GixxbHdXX3X0dH1szVGMShpPiYzAq9A2x9gsLQzLgFfN918FvjP1mVLqiFIq33y/AjgFhFq4XYtUn64jvXkPZcP0+DCOJHRYBDkeaQyr2Gx0FIcWULqRY66RBMakGh1F6yMRoXzYHBKbD9BYVz3g27P0r2KYUuo4gPnrsIs1FpHJgAdQ2G3xY+YupqdExCYz5ORsX2uewvM6W2xOs6IzMQuJ7TzKiZIco6M4pPKKCtLas6iMXGB0FK2fho6/jkNuKZw6MfADSvZaGERkk4gc6uHWr1G3RCQceB24Uyl19ujhL4ExwCQgCPjFRda/V0QyRSSzsrKyP5v+7msd/oh6fIjJ0H2sjibq7BwNNuprdTYFO97HXToJ03MvOJzUKXMZ/5sviRs9dsC31WthUErNV0ql9nD7EDhp/oN/9g//qZ5eQ0T8gY+BXyuldnd77eOqSyvwMjD5IjlWKqUylFIZoaGX3hPV0NxCSv1OioNm6ik8HVB0fBIFLnH4Hd1gdBSH5F24nmoJIjxphtFRtH4SG549aWlX0lrgdvP924EPz28gIh7AGuA1pdS75z13tqgIXccnDlmYp1eHdn5CoNTjmfadwyGagzgZMY/E1hxq9RwN/VJZU0N6815Kh8/Xx9a0i7L0p+NxYIGI5AMLzI8RkQwR+Ye5zY3A5cAdPZyW+oaIZAFZQAjwBwvz9Ko160Na8GCUnsLTYQ2bciMuoij4/C2joziU3C/ew1vaCMrQp2hrF+dmycpKqWpgXg/LM4G7zff/CfzzAuvPtWT7/dXW3sGY2s8p8J9CqpevLTetWdGolEmUvh+BT8E64GGj4zgMj7yPqJWhRI/7zq+spp1jUO1PZu3dSpjUIElLjY6iWUBcXCgLX8jolgPUVVUYHcch1J6uI61pN0dD5yCuFv0/qA0Cg6ow1H/zAe3KlfiZelfa0YVOuRE3MZH/+dtGR3EI2ds/YIi04j9B/+xrvRtUhSHcdJyj/hPw8g8xOopmoYS0aZTLcDzzbT8friNyyV1LHb7E6VO0tT4YVPuUo3/8PrS3GB1DswJxcaFs+AImVbzBmeqT+AeHGR3Jbp1paCCtYQeFIXMZ5+ZhdBzNAQyqPQYA9Ny2TiNoUld3Ut4XujvpYrK3r8VPmhkyXncjaX0z+AqD5jQSx82kQobhcUR3J12MKftD6vEhfrKelEfrG10YNIclLi4cDZtPUtM+6mqrjI5jlxqbmkk58yVFgZfh4qH3lrW+0YVBc2hBGTfgIZ3k6bOTenRox8cESCNeY/WAkVrf6cKgObTECbM5QQgeebo7qSfth9bQhBejpl1tdBTNgejCoDk0cXGhZNg8kpv2cqauxug4dqWlpYXk059TEDADV08fo+NoDkQXBs3hBWTciId0cPjzd4yOYleyd3xEkNTjlq7PRtL6RxcGzeElTpzDKYJwz/3A6Ch2pePAu9TjQ8KMa42OojkYXRg0h+fi6krx8CtIadrD6eoepwQZdJqaGkiu+4IjgbNx9/Q2Oo7mYHRh0JxCyLSb8ZBODm99w+godiH7i/fxk2a8J9xodBTNAenCoDmFkWkzKXeJYMgR3Z0EIFnvUYs/o6fpkYS1/tOFQXMK4uJCReSVpLQe4ER5sdFxDHXmTC0pDTspCJ2Hq5u70XE0B2RRYRCRIBHZKCL55q+BF2jX2W32trXdlseJyFfm9d82TwOqaZdkxGW34CKKom2vGx3FULnb3sFb2vCftMLoKJqDsnSP4RFgs1IqAdhsftyTZqXUOPOt+5U2fwKeMq9fC9xlYR5tEBuRMI5Ct3iCi9f23tiJeeS+zymCScyYb3QUzUFZWhiWAa+a778KXNPXFUVEgLnA6ktZX9N6UhN3NaM78yk5ctDoKIaorT5FStNeSsKvQFxcjY6jOShLC0OYUuo4gPnrsAu08xKRTBHZLSJn//gHA6eVUh3mx+XAiAttSETuNb9GZmVlpYWxNWcVN/tWAI592eM0404vb+sbeEgnIVNvNjqK5sB6nahHRDYBw3t46lf92E60UqpCREYCW0QkCzjTQzt1oRdQSq0EVgJkZGRcsJ02uIWMiOewZxqR5R+jTI8jLoPr/Ar/I+9RLhHEpc0wOormwHr9rVFKzVdKpfZw+xA4KSLhAOavPV5dpJSqMH8tArYB44EqIEBEzhanSEDP7K5ZrDHxGmJUOXkHdhgdxaYqS/NIbsviaNSyQVcQNeuy9KdnLXC7+f7twIfnNxCRQBHxNN8PAWYAOUopBWwFll9sfU3rr8S5t9KmXKnd+ZrRUWyqeOvLAETPvtPgJJqjs7QwPA4sEJF8YIH5MSKSISL/MLdJAjJF5ABdheBxpVSO+blfAD8TkQK6jjmssjCPpuEXGEa23wxGV35Ca2uz0XFsQplMjDj6AVnuY4kaOdroOJqD6/UYw8UopaqBeT0szwTuNt/fCaRdYP0iYLIlGTStJ24TbyHo8y/I3PYeGVfcYnScAVe0fyvxpuOUpvzI6CiaE9AdkZpTSp55LdUE4HLgX0ZHsYnTO1+lSXmSPNf5i6A28HRh0JySq7sHheFLSGvcTeWJcqPjDKj21iYSqzaS5X8ZQwODjI6jOQFdGDSnFT7rB7hLJ0c2vWx0lAF1eNvb+NGE+4TvGx1FcxK6MGhOK2rMJArcEhhW9B5dJ8E5qQNvcYJg0mZeZXQSzUnowqA5tbrE5SSYijm83zmvaThzooTkxq84EnYl7u56JFXNOnRh0Jxa4oI7aVOu1Ox4xegoA6Jo4/MIEDb7XqOjaE5EFwbNqfkFhpHtfzmpletpbKg3Oo5Vqc4ORhS9yzce4xmd1OMZ4Zp2SXRh0Jye9/R7GCqNZH32itFRrKpo94eEqioaU/Upqpp16cKgOb3RUxZT5jKCodmvO9VB6JbdL1GphjJ+vp6QR7MuXRg0pycuLpxIvJmkzjzynOQgdGNVKWPO7CAr9Cr8hvgYHUdzMrowaINC0qL7aFHu1Hyx0ugoVlG44XlcRTFMH3TWBoAuDNqg4BsQSk7QfMbWbOB0bY3RcSyiOtsZXvgO+9zGkZKSbnQczQnpwqANGoGzfsgQaeHQpy8aHcUiRV++zTBTJWfS7qRrhlxNsy5dGLRBI27sLIrdRjIi/5+YOk1Gx7l0X/2dMsKYvFAfdNYGhi4M2uAhwun0e4gzlXLgizVGp7kkVXm7iW8+RG7UTQzx9jQ6juakLJqPwZ60t7dTXl5OS0uL0VGcmpeXF5GRkQ47/ELaFXdS9fUTyO7nYM71Rsfpt5Mbn8ZLeZG8RM+7oA0ciwqDiAQBbwOxQAlwo1Kq9rw2c4Cnui0aA9yklPpARF4BZgF15ufuUErtv5Qs5eXl+Pn5ERsbq/tdB4hSiurqasrLy4mLizM6ziVx8/SmZOTNZBQ9x5GDe0hMd5x5olpqjpFYtZEvhl7FvPDhRsfRnJilXUmPAJuVUgnAZvPjcyiltiqlximlxgFzgSbgs25Nfn72+UstCgAtLS0EBwfrojCARITg4GCH3ysbfdWDNCsPqjc/1XtjO1L4yTO4qk4CZ/+H0VE0J2dpYVgGvGq+/ypwTS/tlwOfKKWaLNxuj3RRGHjO8D32CwwjZ9hSJpzeyPFjpUbH6RNT8xmi8v/JbvfJjB+fYXQczclZWhjClFLHAcxfh/XS/ibgzfOWPSYiB0XkKRG54NE0EblXRDJFJLOystKy1A7C19cXgIqKCpYvX37Rtk8//TRNTf+ut0uWLOH06dMDms+RRS75T9zpoPCjvxgdpU/yP3kWfxromP6QUxRnzb71WhhEZJOIHOrhtqw/GxKRcCAN2NBt8S/pOuYwCQgCfnGh9ZVSK5VSGUqpjNDQ0P5s2q50dnb2e52IiAhWr1590TbnF4b169cTEBDQ720NFmFxqRz0n8W44+9QXXXS6DgXpdqbCT30Ipku6UyftcjoONog0GthUErNV0ql9nD7EDhp/oN/9g//qYu81I3AGqVUe7fXPq66tAIvA45zJLAHJSUljBkzhttvv5309HSWL19OU1MTsbGxPProo8ycOZN3332XwsJCFi1axMSJE7nssss4fPgwAMXFxUybNo1Jkybxm9/85pzXTU1NBboKy8MPP0xaWhrp6ek8++yzPPPMM1RUVDBnzhzmzJkDQGxsLFVVVQA8+eSTpKamkpqaytNPP/3tayYlJXHPPfeQkpLCwoULaW5uBuCZZ54hOTmZ9PR0brrpJpt9/2wteMmv8JVmctY8YXSUiyra+CJBplpOZ/wYN1d9hrk28Cw9XXUtcDvwuPnrhxdpu4KuPYRviUi4Uuq4dO0bXwMcsjAPAL//KJucijPWeKlvJUf489urUnptl5eXx6pVq5gxYwY/+MEP+Nvf/gZ0nea5fft2AObNm8fzzz9PQkICX331FT/60Y/YsmULDz74IPfffz+33XYbzz33XI+vv3LlSoqLi/nmm29wc3OjpqaGoKAgnnzySbZu3UpISMg57fft28fLL7/MV199hVKKKVOmMGvWLAIDA8nPz+fNN9/kxRdf5MYbb+S9997jlltu4fHHH6e4uBhPT0+n7o6KSprMQd+ZpJe/SU31LwkKDul9JVvrbMdv33MclERmzr/O6DTaIGHpvx+PAwtEJB9YYH6MiGSIyD/ONhKRWCAK+Py89d8QkSwgCwgB/mBhHsNFRUUxY8YMAG655ZZvi8H3vvc9ABoaGti5cyc33HAD48aN47777uP48eMA7NixgxUruq5mvfXWW3t8/U2bNvHDH/4QN7eumh4UFHTRPNu3b+faa69lyJAh+Pr6ct111/Hll18CEBcXx7hx4wCYOHEiJSUlAKSnp/P973+ff/7zn99ux1kFLvoVQ6WRg2v+bHSUHpVseoFhnSc4lvYfeHk492eh2Q+LftKUUtXAvB6WZwJ3d3tcAozood1cS7Z/IX35z36gnH9g8OzjIUOGAGAymQgICGD//p7PzO3twKJSql8HHy82/4Cn57+P9bu6un7blfTxxx/zxRdfsHbtWv7nf/6H7Oxspy0QUanTyf5kKuPLXqe66mcEh4QZHelbqq0J/6+eZL+MYfaVejIezXZ0h6WVlZaWsmvXLgDefPNNZs6cec7z/v7+xMXF8e677wJdf7gPHDgAwIwZM3jrrbcAeOONN3p8/YULF/L888/T0dEBQE1N10ihfn5+1Nd/d+rKyy+/nA8++ICmpiYaGxtZs2YNl1122QXzm0wmysrKmDNnDk888QSnT5+moaGhP98Ch+O/9FH8aCLnnd8ZHeUcheufIshUzclJv8Db0zkLs2afdGGwsqSkJF599VXS09Opqanh/vvv/06bN954g1WrVjF27FhSUlL48MOuQzN//etfee6555g0aRJ1dXXfWQ/g7rvvJjo6mvT0dMaOHcu//vUvAO69914WL1787cHnsyZMmMAdd9zB5MmTmTJlCnfffTfjx4+/YP7Ozk5uueUW0tLSGD9+PD/96U+d/uymqKQpHAhexOST71JSmGt0HABMTbWEHfgbu1wmMGdhb5cHaZp1iSNOdZiRkaEyMzPPWZabm0tSUpJBibqUlJSwdOlSDh2yyjF0u2UP32trq60owueFyezzncX0n79ndBzy/vkzRhesYvOs1cybs8DoOJqTEJF9Sqler5DUewyaBgRGjCQn5vtMbdjM/q+2GpqlseIwIwteZbPnfObMmm9oFm1w0oXBimJjY51+b8GZJd3wO2plKF4bfk5rW5sxIZTixNsP0azcCbv+cVxc9FXOmu3pwqBpZl5+gZyY+hvGmPLZ+fb/GpLh+J73ia/bxdbwu0hNTDAkg6bpwqBp3aRccReHfSYwseAZSkqKbLrtzqZaPDb8nHyimH7TL3tfQdMGiC4MmtadCKHfexYvaaPyzfvp6Oj/2FaXqvD1nzC0s5byy/9CaICvzbaraefThUHTzhMck0pe6n8yqXU3296yTZfS8b1rSDy+lk8CVzB7zhU22aamXYguDJrWg7TrHiFvyESm5/+FrAOZva9ggZbqMrzXP8gRoply++N6WG3NcLow2JE77rij1+G1z/f888/z2muvDVCiQczFhRF3vEK7eOD9wQ+orK4emO10tnPiHzfhbmqhdskLDAv0H5jtaFo/6MIwQJRSmEymAd/OD3/4Q2677bYB385g5BsaTe2SvxNnKiX/hVtpaWvvfaX+UIrcl+8ntvkQW0f/N1MmT7fu62vaJXLOAVg+eQROZFn3NYenweLHL9qkpKTk22Epdu3axUMPPcTzzz9Pa2sr8fHxvPzyy/j6+vLoo4/y0Ucf0dzczPTp03nhhRf61H3wyCOPsHbtWtzc3Fi4cCF/+ctf+N3vfoevry8PP/wws2fPZsqUKWzdupXTp0+zatWqi46LpPUudvJVHD72C6YfeJzPnvsRc378PO5urlZ57cOrf09S+btsCLiRxTfpeZw1+6H3GKwsLy+P2267jY0bN7Jq1So2bdrE119/TUZGBk8++SQADzzwAHv37uXQoUM0Nzezbt26Xl+3pqaGNWvWkJ2dzcGDB/n1r3/dY7uOjg727NnD008/ze9//3urvrfBasw1j5Ab9T0W1r3DZ3/7GR2dlu8J5qx9kjHZT/Gl12xm/ehvuOoL2TQ74px7DL38Zz+QYmJimDp1KuvWrSMnJ+fbuRna2tqYNm0aAFu3buWJJ56gqamJmpoaUlJSuOqqqy76uv7+/nh5eXH33Xdz5ZVXsnTp0h7bXXdd12Qu3edX0CwkQtKdz3PkxWauPP4Kn/61gRn3/w0/7wtOUX5hSpH19m9JO/xX9npMJv0//oWXh7v1M2uaBSzaYxCRG0QkW0RMInLBgZlEZJGI5IlIgYg80m15nIh8JSL5IvK2iHhYkscenJ13QSnFggUL2L9/P/v37ycnJ4dVq1bR0tLCj370I1avXk1WVhb33HMPLS0tvb6um5sbe/bs4frrr+eDDz5g0aKe5/49O8eCq6vrt0Nza1bg4kLiPa+QF7OCRWdWc+B/ryb/aFm/XqK16Qz7n7mJtMN/ZYf3HJIfWstQvyEDk1fTLGBpV9Ih4Drgiws1EBFX4DlgMZAMrBCRZPPTfwKeUkolALXAXRbmsRtTp05lx44dFBQUANDU1MSRI0e+LQIhISE0NDT0+SykhoYG6urqWLJkCU8//fQFJ/rRBpCLK6Pv+DvFGb9mSsde/F66nA/ffJ4zzRcfV0kpxYFt73Hqz1NIr9nA5+F3Melnqxni422j4JrWP5bO4JYLvc46NhkoUEoVmdu+BSwTkVxgLnCzud2rwO+Av1uSyV6EhobyyiuvsGLFClpbWwH4wx/+QGJiIvfccw9paWnExsYyadKkPr1efX09y5Yto6WlBaUUTz311EDG1y5EhLilP+f0mFmY3rmHZXm/IPtPL1AQuZywjGUkjEogwMeDtg4T5ceOUrF3HaF5/2JsZy5lEsHBua8ya9Yyo9+Fpl2UVeZjEJFtwMPmKT3Pf245sEgpdbf58a3AFLqKwG6l1Cjz8ijgE6VUam/bs9f5GAYL/b026+yg/POXcd39LOFtRwGoUb7UKx/8pIkg6Zr5rsJlOMfH3EHKVQ/i5e1jZGJtkOvrfAy97jGIyCZgeA9P/Uop9WFfsvSwTF1k+YVy3AvcCxAdHd2HzWraAHN1I3LuPTDnbppLv6Hi4BbaTuQgbQ3Uu/lyPCiW4NT5RIyeSoSLPgFQcxy9FgallKUzhZQDUd0eRwIVQBUQICJuSqmObssvlGMlsBK69hgszGTXrr32WoqLi89Z9qc//YkrrtBj6NglEbxjJhAfM8HoJJpmFbY4XXUvkCAiccAx4CbgZqWUEpGtwHLgLeB2oC97IE5vzZo1RkfQNG0Qs/R01WtFpByYBnwsIhvMyyNEZD2AeW/gAWADkAu8o5TKNr/EL4CfiUgBEAyssiSPI85f7Wj091jTnJ+lZyWtAb7z761SqgJY0u3xemB9D+2K6DpryWJeXl5UV1cTHBysR6ccIEopqqur8fLyMjqKpmkDyGmufI6MjKS8vJzKykqjozg1Ly8vIiMjjY6hadoAcprC4O7uTlxcnNExNE3THJ4+h07TNE07hy4MmqZp2jl0YdA0TdPOYZUhMWxNRCqBo5e4eghdF9cNJvo9Dw76PQ8OlrznGKVUaG+NHLIwWEJEMvsyVogz0e95cNDveXCwxXvWXUmapmnaOXRh0DRN084xGAvDSqMDGEC/58FBv+fBYcDf86A7xqBpmqZd3GDcY9A0TdMuYlAVBhFZJCJ5IlIgIo8YnccaRCRKRLaKSK6IZIvIg+blQSKyUUTyzV8DzctFRJ4xfw8OiojDTiIgIq4i8o2IrDM/jhORr8zv+W0R8TAv9zQ/LjA/H2tk7kslIgEislpEDps/72nO/jmLyE/NP9eHRORNEfFyts9ZRF4SkVMicqjbsn5/riJyu7l9vojcbkmmQVMYRMQVeA5YDCQDK0Qk2dhUVtEB/KdSKgmYCvyH+X09AmxWSiUAm82Poev9J5hv9+LYc2w/SNdQ7mf9CXjK/J5rgbvMy+8Cas3TyD5lbueI/gp8qpQaA4yl67077ecsIiOAnwAZ5il/Xemaz8XZPudXgEXnLevX5yoiQcBv6Zo2eTLw27PF5JIopQbFja45IzZ0e/xL4JdG5xqA9/khsADIA8LNy8KBPPP9F4AV3dp/286RbnTN+LcZmAuso2uq2CrA7fzPm665QKaZ77uZ24nR76Gf79cfKD4/tzN/zsAIoAwIMn9u64ArnPFzBmKBQ5f6uQIrgBe6LT+nXX9vg2aPgX//kJ1Vbl7mNMy7zuOBr4AwpdRxAPPXYeZmzvJ9eBr4LxEOR8oAAAJfSURBVMBkfhwMnFZdE0PBue/r2/dsfr7O3N6RjAQqgZfN3Wf/EJEhOPHnrJQ6BvwFKAWO0/W57cO5P+ez+vu5WvXzHkyFoafZe5zmlCwR8QXeAx5SSp25WNMeljnU90FElgKnlFL7ui/uoanqw3OOwg2YAPxdKTUeaOTf3Qs9cfj3bO4KWQbEARHAELq6Us7nTJ9zby70Hq363gdTYSgHoro9jgQqDMpiVSLiTldReEMp9b558UkRCTc/Hw6cMi93hu/DDOBqESmha77wuXTtQQSIyNk5Rrq/r2/fs/n5oUCNLQNbQTlQrpT6yvx4NV2Fwpk/5/lAsVKqUinVDrwPTMe5P+ez+vu5WvXzHkyFYS+QYD6jwYOug1hrDc5kMemax3QVkKuUerLbU2uBs2cm3E7XsYezy28zn90wFag7u8vqKJRSv1RKRSqlYun6HLcopb4PbAWWm5ud/57Pfi+Wm9s71H+SSqkTQJmIjDYvmgfk4MSfM11dSFNFxMf8c372PTvt59xNfz/XDcBCEQk072ktNC+7NEYfdLHxAZ4lwBGgEPiV0Xms9J5m0rXLeBDYb74toatvdTOQb/4aZG4vdJ2dVQhk0XXGh+Hvw4L3PxtYZ74/EtgDFADvAp7m5V7mxwXm50canfsS3+s4INP8WX8ABDr75wz8/v9v3w6REIaBAIp+xzlpj4bgQhUV3ASDKGZdBxy8d4JkVvxMMqn2aqtu1eXX5lzdO95Qnh0n//WTuVbLe++P6vrNmvx8BmD4p6skAE4QBgAGYQBgEAYABmEAYBAGAAZhAGAQBgCGF+l1LPN3Z0nwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# 1. 定义RNN的参数\n",
    "HIDDEN_SIZE = 30                            # LSTM中隐藏节点的个数。\n",
    "NUM_LAYERS = 2                              # LSTM的层数。\n",
    "TIMESTEPS = 10                              # 循环神经网络的训练序列长度。\n",
    "TRAINING_STEPS = 10000                      # 训练轮数。\n",
    "BATCH_SIZE = 32                             # batch大小。\n",
    "TRAINING_EXAMPLES = 10000                   # 训练数据个数。\n",
    "TESTING_EXAMPLES = 1000                     # 测试数据个数。\n",
    "SAMPLE_GAP = 0.01                           # 采样间隔。\n",
    "\n",
    "\n",
    "# 2. 产生正弦数据函数\n",
    "def generate_data(seq):\n",
    "    X = []\n",
    "    y = []\n",
    "    # 序列的第i项和后面的TIMESTEPS-1项合在一起作为输入；第i + TIMESTEPS项作为输出。\n",
    "    # 即用sin函数前面的TIMESTEPS个点的信息，预测第i + TIMESTEPS个点的函数值。\n",
    "    for i in range(len(seq) - TIMESTEPS):\n",
    "        X.append([seq[i: i + TIMESTEPS]])\n",
    "        y.append([seq[i + TIMESTEPS]])\n",
    "    return np.array(X, dtype=np.float32), np.array(y, dtype=np.float32)  \n",
    "\n",
    "\n",
    "# 3. 定义网络结构和优化步骤\n",
    "def lstm_model(X, y, is_training):\n",
    "    # 使用多层的LSTM结构。\n",
    "    cell = tf.nn.rnn_cell.MultiRNNCell([\n",
    "        tf.nn.rnn_cell.LSTMCell(HIDDEN_SIZE)\n",
    "        for _ in range(NUM_LAYERS)]) \n",
    "\n",
    "    # 使用TensorFlow接口将多层的LSTM结构连接成RNN网络并计算其前向传播结果。\n",
    "    outputs, _ = tf.nn.dynamic_rnn(cell, X, dtype=tf.float32)\n",
    "    # outputs是顶层LSTM在每一步的输出结果，它的维度是[batch_size, time ,\n",
    "    # HIDDEN_SIZE]。在本问题中只关注最后一个时刻的输出结果。\n",
    "    output = outputs[:, -1, :]\n",
    "\n",
    "    # 对LSTM网络的输出再做加一层全链接层并计算损失。注意这里默认的损失为平均\n",
    "    # 平方差损失函数。\n",
    "    predictions = tf.contrib.layers.fully_connected(\n",
    "        output, 1, activation_fn=None)\n",
    "    \n",
    "    # 只在训练时计算损失函数和优化步骤。测试时直接返回预测结果。\n",
    "    if not is_training:\n",
    "        return predictions, None, None\n",
    "        \n",
    "    # 计算损失函数。\n",
    "    loss = tf.losses.mean_squared_error(labels=y, predictions=predictions)\n",
    "\n",
    "    # 创建模型优化器并得到优化步骤。\n",
    "    train_op = tf.contrib.layers.optimize_loss(\n",
    "        loss, tf.train.get_global_step(),\n",
    "        optimizer=\"Adagrad\", learning_rate=0.1)\n",
    "    \n",
    "    return predictions, loss, train_op\n",
    "\n",
    "\n",
    "# 4. 定义训练方法\n",
    "def train(sess, train_X, train_Y):\n",
    "    # 将训练数据以数据集的方式提供给计算图\n",
    "    ds = tf.data.Dataset.from_tensor_slices((train_X, train_Y))\n",
    "    ds = ds.repeat().shuffle(1000).batch(BATCH_SIZE)\n",
    "    X, y = ds.make_one_shot_iterator().get_next()\n",
    "    \n",
    "    # 定义模型，得到预测结果、损失函数，和训练操作。\n",
    "    with tf.variable_scope(\"model\"):\n",
    "        _, loss, train_op = lstm_model(X, y, True)\n",
    "        \n",
    "    sess.run(tf.global_variables_initializer())\n",
    "    for i in range(TRAINING_STEPS):\n",
    "        _, l = sess.run([train_op, loss])\n",
    "        if i % 1000 == 0:\n",
    "            print(\"train step: \" + str(i) + \", loss: \", str(l))\n",
    "            \n",
    "\n",
    "# 5. 定义测试方法\n",
    "def run_eval(sess, test_X, test_y):\n",
    "    # 将测试数据以数据集的方式提供给计算图。\n",
    "    ds = tf.data.Dataset.from_tensor_slices((test_X, test_y))\n",
    "    ds = ds.batch(1)\n",
    "    X, y = ds.make_one_shot_iterator().get_next()\n",
    "    \n",
    "    # 调用模型得到计算结果。这里不需要输入真实的y值。\n",
    "    with tf.variable_scope(\"model\", reuse=True):\n",
    "        prediction, _, _ = lstm_model(X, [0.0], False)\n",
    "    \n",
    "    # 将预测结果存入一个数组。\n",
    "    predictions = []\n",
    "    labels = []\n",
    "    for i in range(TESTING_EXAMPLES):\n",
    "        p, l = sess.run([prediction, y])\n",
    "        predictions.append(p)\n",
    "        labels.append(l)\n",
    "\n",
    "    # 计算rmse作为评价指标。\n",
    "    predictions = np.array(predictions).squeeze()\n",
    "    labels = np.array(labels).squeeze()\n",
    "    rmse = np.sqrt(((predictions - labels) ** 2).mean(axis=0))\n",
    "    print(\"Root Mean Square Error is: %f\" % rmse)\n",
    "    \n",
    "    # 对预测的sin函数曲线进行绘图。\n",
    "    plt.figure()\n",
    "    plt.plot(predictions, label='predictions')\n",
    "    plt.plot(labels, label='real_sin')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    \n",
    "# 6. 生成数据并训练、验证\n",
    "# 用正弦函数生成训练和测试数据集合。\n",
    "# numpy.linspace函数可以创建一个等差序列的数组，它常用的参数有三个参数，\n",
    "# 第一个参数表示起始值，第二个参数表示终止值，第三个参数表示数列的长度。\n",
    "# 例如linespace(1, 10, 10)产生的数组是arrray([1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) \n",
    "test_start = (TRAINING_EXAMPLES + TIMESTEPS) * SAMPLE_GAP\n",
    "test_end = test_start + (TESTING_EXAMPLES + TIMESTEPS) * SAMPLE_GAP\n",
    "train_X, train_y = generate_data(np.sin(np.linspace(\n",
    "    0, test_start, TRAINING_EXAMPLES + TIMESTEPS, dtype=np.float32)))\n",
    "test_X, test_y = generate_data(np.sin(np.linspace(\n",
    "    test_start, test_end, TESTING_EXAMPLES + TIMESTEPS, dtype=np.float32)))\n",
    "\n",
    "with tf.Session() as sess:\n",
    "    train(sess, train_X, train_y)\n",
    "    run_eval(sess, test_X, test_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "可以看到，预测得到的结果和真实的$sin$函数几乎是重合的，也就是说通过RNN可以非常好地预测$sin$函数的取值。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
